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Abstract — Subsurface radioactivity may be due to transport 
of radionuclides from a contaminated surface into the solid 
volume, as occurs for radioactive fallout deposited on soil, or from 
fast neutron activation of a solid volume, as occurs in concrete 
blocks used for radiation shielding. For purposes including 
fate and transport studies of radionuclides in the environment, 
decommissioning and decontamination of radiation facilities, and 
nuclear forensics, an in situ, nondestructive method for ascertain- 
ing the subsurface distribution of radioactivity is desired. The 
method developed here obtains a polynomial expression for the 
radioactivity depth profile, using a small set of gamma-ray count 
rates measured by a collimated detector directed towards the 
surface at a variety of angles with respect to the surface normal. 
To demonstrate its capabilities, this polynomial method is applied 
to the simple case where the radioactivity is maximal at the 
surface and decreases exponentially with depth below the surface, 
and to the more difficult case where the maximal radioactivity 
is below the surface. 

Index Terms — Gamma ray detection, nuclear measurements, 
radiation imaging, radioactive materials 



I. Introduction 

Nondestructive methods to determine the extent of sub- 
surface radioactivity rely on the attenuation of characteris- 
tic gamma-rays emitted from the subsurface. Assuming the 
radionuclide distribution is reasonably uniform in the two 
dimensions parallel to the surface, the problem reduces to 
finding the (one-dimensional) depth profile of the radioactivity. 
To accomplish this, two general approaches have been taken: 
(i) A depth distribution function is assumed (from physical 
considerations), and spectroscopic measurements are made to 
determine the best-fit values for the function parameters, (ii) 
No depth distribution is assumed; rather, the subsurface vol- 
ume is divided into voxels, and spectroscopic measurements 
are made to determine the activity of each voxel. 

Often there is good reason to assume a particular depth 
distribution function. For example, Smith and Elder | con- 
sider several parameterized functions based on physical pro- 
cesses related to radionuclide mobility in soils. The review by 
MacDonald et al. ||2l discusses three spectroscopic methods 
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used to obtain the function parameters: (i) The "multiple peak 
method", which uses the fact that the scatter cross-section of a 
gamma-ray is a function of its energy (so for example, the two 
parameters of an exponential decay function representing the 
radionuclide distribution are calculated from measurements of 
the count rates of two characteristic gamma-rays of different 
energy), (ii) The "peak to valley method", which uses the 
fact that the ratio of counts in a characteristic gamma-ray 
peak to counts in the adjacent "valley" of slightly lower 
energy (due to forward scattering of the characteristic gamma- 
rays) decreases with increased attenuation of the characteristic 
gamma-rays (e.g. Gering et al. f3l who use a two-parameter 
Lorentz function for ^"^^Cs in soil), (iii) The "lead plate 
method", where an absorbing lead plate inserted between the 
detector and the material surface allows the flux contribution 
from large angle (with respect to the surface normal) to be 
distinguished from that from small angle (e.g., Korun et al. 
|4| who use a three-parameter Gaussian function for ^^^^Cs and 
^•^^Cs in soil; Dewey et al. [5J who fit a variety of analytical 
functions to simulated measurements for ^'^^Cs in soil; Benke 
and Kearfott |6|,|7| who consider planar radiation sources at 
various depths). 

When a particular depth distribution cannot be assumed, or 
none of a set of parameterized functions can be well matched 
to the spectroscopic measurements, the subsurface volume 
may be (virtually) divided into voxels, and the activity of 
each voxel determined by some variation of the following 
procedure. Count rates of a characteristic gamma-ray are 
obtained by a collimated detector directed at (or otherwise 
sensitive only to gamma-rays emitted at) particular angles 6 
with respect to the surface normal. By obtaining count rates 
for a number of different detector orientations, the response 
matrix equation C = RA can be constructed, where element 
Ci of vector C is the count rate obtained for angle Oi, the 
detector response matrix element is the probability that a 
gamma-ray emitted from voxel j will contribute to the count 
rate Ci (so Vij accounts for scattering of that gamma-ray en 
route to the detector), and element aj of vector A is the activity 
of voxel j. Given the vector C, the matrix equation is solved 
for the activity vector A. This "voxel approach" is taken by, 
e.g., Benke et al. H,||9|; Whetstone et al. ITOl; Charles et al. 

frn. 

These approaches to finding the radioactivity depth profile 
are not entirely satisfactory for the following reasons. In the 
case of the "parameterized function approach", the assumption 
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of a particular depth distribution function is risky in general, 
even when best-fit parameter values are found that seem to 
bear out the assumption. By the nature of the measurements, 
the candidate functions are limited to those with very few 
parameters, so reasonable best-fit parameter values may be 
found for more than one of those functions. In the case of 
the "voxel approach", the limited number of spectroscopic 
measurements greatly limits the number of voxels into which 
the subsurface is divided. In any event, the calculated activity 
Uj (gamma-ray emissions per unit time) is evenly distributed 
over the entire volume of voxel j, so giving the calculated 
distribution a spatial resolution no better than the voxel size. 

The new approach presented below is motivated by the 
expectation that the "true" depth profile function is smooth 
and continuous. Thus it can be closely approximated by a 
polynomial expression, with more or fewer terms depending 
on the complexity of the "true" function (for example, more 
terms are needed to approximate a Gaussian-like function than 
an exponential-like function). The success of this method is 
due to the novel way the coefficients of the terms in the 
polynomial are obtained from the spectroscopic measurements. 

II. Calculation of the polynomial profile 

FUNCTION 

For convenience, the surface is at z = and the radioactive 
material lies in the region z > 0. In what follows, p{z)dV 
is the activity (number of gamma-rays emitted per second) of 
the infinitesimal volume dV located in the subsurface layer 
[z,z + d2]; thus the distribution p(z) is the activity per unit 
volume at depth z. 

This "polynomial method" is most easily developed by 
considering the idealized case of a highly collimated, one- 
pixel gamma-ray detector pointing at the material surface. 
Due to the collimation, only those characteristic gamma-rays 
emitted from a thin "pencil" of material will be counted at the 
detector. Unless all radioactivity is on the surface, the detector 
counts will provide an underestimate of the total activity, 
due to attenuation (scattering) of the subsurface emissions. 
Different orientations of the detector with respect to the surface 
normal will thus produce different count rates, due to the 
different orientations of the pencil of material with respect 
to the subsurface radioactivity profile. As shown below, that 
radioactivity profile can be ascertained from the variation in 
count rate with detector orientation. 

Consider that the highly collimated, one-pixel detector is 
pointing at the surface at an angle 9 from the surface nor- 
mal. The collimation ensures that all recorded characteristic 
gamma-rays originate from the thin pencil of material angled 
by 9 with respect to the surface normal. Note that a gamma-ray 
emitted at depth z and recorded at the detector has travelled 
a distance z/ cos 6* to the surface. Thus the count rate d{9) 
recorded at the detector is given by 
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where the (unknown) function p(z) is the activity per unit 
volume at depth z; /i is the attenuation coefficient for the char- 
acteristic gamma-ray in the homogeneous material (attenuation 
in air is ignored); R is the normal distance of the pixel center 
above the surface; A is the cross-sectional area of the pencil 
volume, as well as the area of the detector pixel; and e is the 
detector efficiency (possibly multiplied by factors related to the 
system geometry and instrument characteristics). Note that the 
integration is over all infinitesimal volumes dV = dz- Aj cos 9 
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comprising the pencil volume. The ratio A 
is the fraction of emissions from the volume dV that would 
reach the detector in the absence of attenuation, while the 
exponential factor accounts for the attenuation. 

When count rates are obtained for several different detector 
orientations (several different 9 values), a set of these integral 
equations is created, which in principle can be solved for the 
unknown radioactivity profile p{z). The approach taken here 
is to approximate p{z) by the polynomial p*{z) of degree N, 

N 

P*{z) — ^ Cj ■ Z^ — Cq + CiZ + C2Z^ + . . . + CnZ^ (2) 

i=o 

where the coefficients Cj are determined from the detector 
count rates in the following manner. 

Each detector orientation with respect to the surface normal 
produces a linear equation 
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where di is the count rate d{9i), as the index i distinguishes 
the different detector orientations; and the coefficients 
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(4) 

Note that the limit of integration ( is the depth below which 
no radioactivity is present; this is a parameter that must be 
set (i.e., assumed) prior to the calculations. In the case of n 
detector orientations, n equations like Eq. ^ are created. Two 
additional linear equations are needed to force the function 
p*{z) to go smoothly to zero at z ~ (. These are 
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where the coefficients K„j — [this is just the equation 
p*(C)=0], and 
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3=0 
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where the coefficients Kn+ij ~ j ■ C"*"^ [this is the equation 
dp*{z)/dz = at z = which ensures the slope of p*{z) is 
zero 2X z — Q. 

The set of n + 2 linear equations is then solved for the set 
{cj}, thus producing the radioactivity profile p*{z) according 
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to Eq. (|2]i. This function reproduces the set {di} of count-rate 
measurements, and goes smoothly to zero at z ~ 

Note that the total activity At (number of characteristic 
gamma-rays emitted per second) of a "core" of cross-sectional 
area A (taken perpendicular to the surface) is 
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The activity of any slice of that core is obtained by integrating 
p*{z) over the thickness of that slice. 

A computational caution: When the integrals in Eq. (|4]i 
are evaluated to obtain the set {Kij}, it may be found that 
those values range over several orders of magnitude. This 
makes it exceedingly difficult to solve the system of n + 2 
equations for the set {cj} by conventional methods. Thus an 
analytical, "term elimination sequence" method for solving a 
square system of linear equations was derived for this purpose. 
It is presented in the appendix in a form that is very easy to 
write into a computer code. All of the polynomial functions 
plotted in the following section were obtained by use of this 
method incorporated into a simple Mathematica code. 

After this exposition of the polynomial method, it is inter- 
esting to briefly consider the "voxel approach" to determining 
the radioactivity depth profile mentioned in Section I. The set 
{di} of count-rate measurements is related to the set {oj} of 
voxel activities by a set of linear equations of the sort 
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where the + 1 voxels are each of size A 
dimensionless coefficients 
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Clearly the number of voxels cannot exceed the number of 
measurements (detector orientations), and the voxel thickness 
5z ^ dz must be sufficient that the array of voxels fills 
the assumed radioactive volume. Of course, the polynomial 
method completely avoids discretization of that volume. 

III. Application of the polynomial method 

Unfortunately, the functional requirements on p*{z) that 
it be continuous, go smoothly to zero at some depth C, 
and reproduce the spectroscopic measurements {di}, do not 
automatically ensure that the unknown radioactivity profile 
function p{z) is reproduced. The finite number of terms in 
the polynomial (limited by the number of measurements taken) 
may be insufficient to produce a complex function, or the finite 
depth ( at which radioactivity is presumed negligible may be 
poorly chosen. Indeed, the exponential decay function, which 
is a typical radioactivity profile, is actually a polynomial with 
an infinite number of terms that goes to zero at infinity. 
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Fig. L Fourth-degree polynomial representations of the radioactivity depth 
profile p (z) = exp [—2/0.003] (the thickest curve). The odd curve is p* (z) 
that terminates at ( = 0.01; the p* (z) curves for C = 0.015, 0.02, 0.025 
track the p (z) curve fairly well. All the p* (z) curves are calculated from 
the count-rate measurements {d. (i^)}. 



But fortunately, non-physical (nonsensical) p*{z) functions 
are easy to distinguish visually. In addition to the functional 
requirements on p*{z), the areas under the curves p*{z) and 
p{z) are approximately equal [see Eq. ^]. Thus either the 
curve p* (z) closely tracks the curve p{z), or the constraints are 
accommodated by anomalous undulations of the curve p*{z). 
Examples of this behavior are evident in Fig. [T] where the ra- 
dioactivity profile p{z) = cxp[— z/0.003] is approximated by 
a 4th degree polynomial calculated from count-rate measure- 
ments at three detector orientations 9 — 0, 7r/16, 7r/8. The four 
curves are p*{z) obtained for ( = 0.01,0.015,0.02,0.025, 
respectively. Clearly the curve obtained for ( = 0.01 is 
physically unlikely due to its hump, while the curve obtained 
for ( = 0.025 actually goes slightly negative. In contrast, 
the curves obtained for ( = 0.015 and 0.02 look physically 
reasonable. Figures |2] and [3] show the 5th degree polynomial 
p*{z) for C, ~ 0.015 and 0.02, respectively. Each figure 
shows three curves (impossible to distinguish since they're 
essentially identical); each curve is calculated from count- 
rate measurements at detector orientations 9 — 0, Lp,2(p,3(p, 
where (p is 7r/20, n/lQ, 7r/12 for the three curves, respectively. 
The p*{z) curves in Fig. |2] are physically unlikely (due to 
the hump), while there is nothing objectionable about the 
p*{z) curves in Fig. [3j Thus any one of the three 5th degree 
polynomials (for C = 0.02) shown in Fig. [3] should be a good 
match to the unknown profile function p{z), as in fact is the 
case. 

Note that the sets {di} [AKA {d{9i)}] of count-rate mea- 
surements used to produce the polynomial curves in these (and 
following) figures are obtained by solving Eq. ([T]i for various 
values of the detector orientation 9. To specify those angles, it 
is convenient to use the notation |rf (^y^) } j_Q 
ample) or just jd (ij^) } when the degree N of the polynomial 
is otherwise noted (the degree of the polynomial is the number 
of detector orientations plus 1). Also, the calculations use the 
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Fig. 2. Fifth-degree polynomial representations of the radioactivity depth 
profile p{z) = exp [—2/0.003] (the thickest curve). The p* (z) curve is 
actually three indistinguishable curves (for f = 0.015) calculated from count- 
rate measurements {d(i^)},{d(iYjr)},{d(i^)}, respectively. 



Fig. 4. Fourth-degree polynomial representations of the radioactivity depth 
profile p(z) = 2 exp [—2/0.003] (the thickest curve). The p* (2) curve is 
actually three indistinguishable curves (for = 0.025) calculated from count- 
rate measurements {d. (i-^)} , {d (^ij^)^ , {d (i^)}, respectively. 




Fig. 3. Similar to Fig. |2] but the three indistinguishable p* (2) curves 
terminate at = 0.02 rather than at = 0.015. The absence of the non- 
physical hump evident in the curves of Fig. |2] suggests that the 5th degree 
polynomial p* (2) for = 0.02 is a good representation of the radioactivity 
depth profile p (2). 



attenuation coefficient /i = 18 appropriate for the 662 
keV gamma-ray emitted by ^'^^Cs in concrete; detector-surface 
separation R = 0.02 m; detector pixel area A — (0.005)^ m^; 
detector efficiency e = 1. 

A more challenging case is that of the radioactivity profile 
p{z) = z exp [— z/0.003], which goes steeply to zero at the 
material surface. A profile of this general shape may occur for 
neutron activated concrete shielding, due to thermalization and 
subsequent capture of fast neutrons incident on the concrete 
blocks (Kimura et al. |12|, Masumoto et al. |fT3l . Wang et 
al. HU). Figure H] shows three (impossible to distinguish) 
4th degree polynomials p*{z) for ( = 0.025 (other, sim- 
ilar calculations showed the appropriate value for to be 
0.02 < C, < 0.03) calculated from count-rate measurements 



Fig. 5. Similar to Fig. |4] but the three p* (2) curves are fifth-degree 
polynomials. 



{d , {d [ijg)} , {d (ij^)}, respectively. These results 

show convincingly that the radioactivity is maximal below the 
surface, but one might wonder about the behavior very close 
to the surface, as a 4th degree polynomial has too few terms 
to effect a tight bend down to a very low surface value. Figure 
|5] shows the corresponding 5th degree polynomials p*{z), all 
of which go nearly to zero at the material surface. 

It is interesting to consider higher-degree (beyond 5th) 
polynomial approximations to the profile function p{z) — 
z exp [— z/0.003] for their cautionary tales, even though it 
may be difficult to get more than four distinct measurements 
di in practice. Figure |6] shows the corresponding 6th degree 
polynomials p*{z). That calculated from the measurement set 
{d[i^)} tracks the p{z) curve very nicely; the others are 
physically impossible! Figures|2]and|8]show the corresponding 
7th and 8th degree polynomials p*{z), respectively. Clearly, 
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Fig. 6. Similar to Figs. |4] and [5] but the three p* {z) curves are sixth- 
degree polynomials. The wildly undulating curves calculated from measure- 
ments {<j'(*^)} and {(^(jj^)} are clearly non-physical. In contrast, the 
reasonable behavior of the p* (z) curve calculated from the measurements 
{d (i^) } suggests that it is a good representation of the radioactivity depth 
profile p (2). 
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Fig. 7. Similar to Figs.|4](6] but the three p* (z) curves are .seveni/i-degree 
polynomials. The wildly undulating curve is calculated from the measurements 
{d (note that it closely resembles the corresponding curve in Fig. |6). 

use of the set {d (ij^)} is preferable, presumably due to its 
wider range of values (of course, all detector orientations must 
have Q < 7r/2). 

The polynomial method also handles the case of a uniform 
radioactivity distribution p{z) — 1, as occurs for natural ra- 
dioactivity in soils for example. Figure |9] shows the 4th degree 
polynomials p*{z) for ( — 0.15,0.5,1,2 calculated from 
count-rate measurements {c?(*y^)}- From visual inspection 
of the curves, it is clear that the polynomial p*{z) —?■ p{z) as 

C — 00. 

IV. Use of multi-pixel detectors 

While a "highly collimated, one-pixel gamma-ray detector" 
is used in Section II to simplify the presentation of the method. 




0.000 0.005 0.010 0.015 0.020 0.025 



z 

Fig. 8. Similar to Figs. |4}0 but the two p* (z) curves are eighth- 
degree polynomials. Again, the wildly undulating curve is calculated from 
the measurements {d. {i-^)}, while the reasonable curve is calculated from 
the measurements { d ( j ) ]• . 
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Fig. 9. Fourth-degree polynomial representations of the uniform radioactivity 
depth profile p (2) = 1 (the thickest curve). The p* (2) curves for (J = 
0.15, 0.5, 1, 2 are calculated from the count-rate measurements {d {ijg)}- 



in practice a collimated detector with many pixels is needed 
to accumulate sufficient counts in a reasonable time. Recall 
that collimation is necessary so that each pixel records only 
characteristic gamma-rays emitted at a particular angle 6 with 
respect to the surface normal. This is perhaps most easily 
accomplished by a collimator that is a thick sheet of lead with 
a set of parallel holes (say, M in number) drilled through it. 
This effectively converts a standard detector into an array of 
M collimated, single-pixel detectors (the single pixel is the 
size of the hole cross-section). Of course, there being only 
a single detector behind the perforated lead collimator, only 
a total gamma-ray count is acquired; that is, the sum of the 
contributions of the M single-pixel detectors is obtained. 

Use of such a detector is easily accommodated by the 
mathematical development in Section II. In particular, the 
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count rate di [AKA d{6i), to make clear that the subscript 
i identifies the particular value of 9] in Eq. Q is then the 
total counts recorded at the detector (that is, the sum of the 
contributions of the M single-pixel detectors) over the time 
period, and Eq. ^ for the coefficients Kij becomes 
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(10) 



where i?^^™^ is the normal distance of pixel m above the 
surface. 

V. Incorporation of known surface conditions 

The polynomial method is intended for use when the "true" 
radioactivity depth profile p (z) is unknown. However, it is 
conceivable that the activity at the material surface is known 
(via a swipe test or a shallow core, for example), thus giving 
a value for the coefficient cq. Then Eq. Q is written 
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which is conveniently expressed as 



N 

i=i 



(11) 



(12) 



where 6i = di — KiQCo and the index i runs from 1 to N, 
so that the term elimination sequence method (see appendix) 
can be used to solve the system of equations for the co- 
efficients ci, . . . , cjv- In the case of the radioactivity profile 
p{z) — z exp [— z/0.003] considered above, prior knowledge 
that Co = would allow Eq. (fTZt to be used, with count-rate 
measurements {d ((i — 1) Y^)}j_;^ ^_2' say. 

Similarly, the slope ci of the radioactivity profile at the sur- 
face may be known as well (for example, buried radionuclides 
would have co = ci = 0), so Eq. (|3]l is expressed as 



N 



(13) 



where 6i = di — KioCQ — KuCi and the index i runs from 
2 to N, and again the term elimination sequence method of 
solution can be used. 

VI. Concluding remarks 

The expectation that the radioactivity depth profile is con- 
tinuous and smooth leads to the "polynomial method" for 
determining that profile. Values for the coefficients of the 
terms in the polynomial are obtained in a clever way from 
a set of characteristic gamma-ray count-rate measurements 
such that the polynomial function exactly reproduces the 
measurements. Practical virtues of this approach include (i) 
it is non-destructive, (ii) no a priori assumptions regarding 



the functional form of the radioactivity depth profile are 
needed, and (iii) very few spectroscopic measurements are 
needed. In contrast, current non-destructive approaches assume 
a particular radioactivity depth profile function (which is then 
fitted to the measurements) or divide the subsurface volume 
into a small number of voxels (which cannot exceed the 
number of distinct measurements). 

The polynomial method is particularly well suited to as- 
certaining the depth distribution of subsurface radioactivity 
due to neutron activation of radiation shielding (Kimura et 
al. |FT2l, Masumoto et al. |13|, Wang et al. |14|, Zagar 
and Ravnik 115|) and structural materials (Nakanishi et al. 
lfT6l ). Thus important applications are to decommissioning and 
decontamination of accelerators and other radiation facilities 
(e.g., enabling separation of radioactive portions from larger 
volumes), and nuclear forensics (e.g., determining the neutron 
fluence and energy spectrum after a nuclear event). 

Appendix: Term elimination sequence method for 
solution of systems of linear equations 

After obtaining the set {di} of count-rate measurements 
and calculating the set {Kij} of coefficients, the system of 
equations of the sort Eq. (O must be solved for the set {cj}. 
Note that 
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This last equation can be written 
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Clearly this two-step procedure can be repeated until the 
number of terms in the summation is reduced to one. This 
suggests the following method of solution. 

Consider the (square) system of linear equations of the sort 
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where the superscripts on the quantities Qi and Qij 
are intended only to distinguish the sets l^^^^ Q^j^'' |, 

l^i^ ^KQij' |, etc., as needed below. The quantities g. 



(TV) 



and Q'^'' correspond to di and Kij, respectively, so the set 
I a,- , Qif''' r" is already known. Then according to 
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the two-step procedure above. 
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this way all 

L '-^ J 0<iJ<Af-l L -'J ij=0 

calculated. Then cq, ci, . . . , cjv, respectively, are obtained 
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